Spatial Prediction with INLA

SPDE Toy Example

Toy dataset from Blangiardo and Cameletti (2015).

# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy,
     pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toymesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]),
                        max.edge = c(0.1, 0.2))
plot(toymesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5),
       pch = 19)

# Initialize SPDE projector matrix.
toyA <- inla.spde.make.A(toymesh, as.matrix(SPDEtoy[,c('s1', 's2')]))

# Initialize exponential covariance structure for SPDE.
toyspde <- inla.spde2.matern(toymesh, alpha = 2)

toy_formula <-  y ~ -1 + intercept + f(spatial_field, model = toyspde)
toy_data <- list(
    y = SPDEtoy$y,
    intercept = rep(1, nrow(SPDEtoy)),
    spatial_field = seq_len(nrow(SPDEtoy))
  )

# Fit the model with INLA.
# Can't find the intercept??
toyfit <- inla(
  toy_formula,
#  y ~ f(spatial_field, model = toyspde),
  data = toy_data,
#  control.predictor = list(A = toyA, compute = TRUE)
#  control.predictor = list(compute = TRUE)
)

Bei Dataset

Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).

# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

bei_corners <- vertices.owin(Window(bei))
bei_domain <- cbind(bei_corners$x, bei_corners$y)
bei_full_mesh <- inla.mesh.2d(cbind(bei$x, bei$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_full_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei, pch = '.', col = 'blue')

bei_full_spdf <- as.SpatialPoints.ppp(bei)
# CHECK PRIORS!
matern_full <- inla.spde2.pcmatern(bei_full_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = matern_full) + Intercept
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf)

lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
n_quads <- 10
botleft <- cbind(runif(n_quads, 0, 950), runif(n_quads, 0, 450))
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
    cbind(
      botleft[r, 1] + c(0, 0, 50, 50),
      botleft[r, 2] + c(0, 50, 50, 0)
    )
  )})
bei_win <- do.call(
  union.owin,
  apply(botleft, 1, function(x){return(
    owin(x[1] + c(0, 50), x[2] + c(0, 50))
  )})
)
bei_hole <- bei[complement.owin(bei_win)]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)

plot(bei_hole, main = 'Observed Subregion', pch = '.', cols = 'black')

bei_hole_mesh <- inla.mesh.2d(cbind(bei_hole$x, bei_hole$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_hole_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei_hole, pch = '.', col = 'blue')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
# CHECK PRIORS!
matern_hole <- inla.spde2.pcmatern(bei_hole_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = matern_hole) + Intercept
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf)

lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

bei_samp_mesh <- inla.mesh.2d(cbind(bei_samp$x, bei_samp$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_samp_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
plot(Window(bei_samp), border = 'blue', add = TRUE)
points(bei_samp, pch = '.', col = 'blue')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
# CHECK PRIORS!
matern_samp <- inla.spde2.pcmatern(bei_samp_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = matern_samp) + Intercept
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf)

lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

References

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.

Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.